act_site = [1 2 4 7 8];
pas_site = [3 5 6 9];
con_site = 10;

for site = 1:9
    for layer = 1:4
        for time = 1:4
            cur_fb_act = proj_meta(site).rd(layer,time).act(:,1:7500);
            cur_pb_act = proj_meta(site).rd(layer,time).act(:,7501:15000);

            % generate binary activity during feedback.
            cur_fb_bi = proj_meta(site).rd(layer,time).velM_smoothed(1:7500);
            cur_fb_bi = -cur_fb_bi; % make the velocity of forward-running positive.
            cur_fb_bi(cur_fb_bi > 0.001) = 1;
            cur_fb_bi(cur_fb_bi <= 0.001) = 0;
            
            % generate binary activity during playback.
            cur_pb_bi = proj_meta(site).rd(layer,time).velM_smoothed(7501:15000);
            cur_pb_bi = -cur_pb_bi; % make the velocity of forward-running positive.
            cur_pb_bi(cur_pb_bi > 0.001) = 1;
            cur_pb_bi(cur_pb_bi <= 0.001) = 0;
            
            % generate the frame id of ps onset.
            cur_ps_ons = find(diff(proj_meta(site).rd(layer,time).ps_id(1:7500))>1)+1;
            % remove early & late ps.
            cur_ps_ons(cur_ps_ons > 7425) = [];
            cur_ps_ons(cur_ps_ons < 6) = [];
            
            % generate the frame id of playback halt onset.
            cur_pb_ons = find(diff(proj_meta(site).rd(layer,time).ps_id(7501:15000))>1)+1;
            % remove early & late pb.
            cur_pb_ons(cur_pb_ons > 7425) = [];
            cur_pb_ons(cur_pb_ons < 8) = [];
            
            rps_nbr = 0; % number of real ps
            rps_id = [];
            for ps_id = 1:length(cur_ps_ons)
                if sum(cur_fb_bi(cur_ps_ons(ps_id)-5:cur_ps_ons(ps_id)+15)) == 21
                    rps_nbr = rps_nbr+1;
                    rps_id(rps_nbr) = ps_id;
                    all_rps_nbr(site,time) = rps_nbr;
                end
            end
            
            rpb_nbr = 0; % number of real pb halt
            rpb_id = [];
            for pb_id = rps_id % rpb is selected based on rps
                if sum(cur_pb_bi(cur_pb_ons(pb_id)-7:cur_pb_ons(pb_id)+22)) == 0
                    rpb_nbr = rpb_nbr+1;
                    rpb_id(rpb_nbr) = pb_id;
                    all_rpb_nbr(site,time) = rpb_nbr;
                end
            end
            
            if rps_nbr >= 10
                rps_id = rps_id(1:10);
%             else
%                 rps_id = [];
            end
            
            if rpb_nbr >= 10
                rpb_id = rpb_id(1:10);
%             else
%                 rps_id = [];
            end
            
            for cell_id = 1:size(cur_fb_act,1)
                rps_resp_cnt = 0;
                rps_resp = [];
                for rps_act_nbr = 1:length(rps_id)
                    rps_frame_id = cur_ps_ons(rps_id(rps_act_nbr));
                    rps_resp_cnt = rps_resp_cnt+1;
                    rps_resp(rps_resp_cnt,:) = cur_fb_act(cell_id,rps_frame_id-7:rps_frame_id+75);
%                     prps_resp(rps_resp_cnt) = cur_fb_act(cell_id,rps_frame_id+8:rps_frame_id+22);
                end
                rps_resp_by_cell{site,layer,time,cell_id} = rps_resp;
%                 prps_resp_by_cell{site,layer,time,cell_id} = prps_resp;
            end
            
            for cell_id = 1:size(cur_pb_act,1)
                rpb_resp_cnt = 0;
                rpb_resp = [];
                for rpb_act_nbr = 1:length(rpb_id)
                    rpb_frame_id = cur_pb_ons(rpb_id(rpb_act_nbr));
                    rpb_resp_cnt = rpb_resp_cnt+1;
                    rpb_resp(rpb_resp_cnt,:) = cur_pb_act(cell_id,rpb_frame_id-7:rpb_frame_id+75);
%                     prpb_resp(rpb_resp_cnt) = mean(cur_pb_act(cell_id,rpb_frame_id+8:rpb_frame_id+22));
                end
                rpb_resp_by_cell{site,layer,time,cell_id} = rpb_resp;
%                 prpb_resp_by_cell{site,layer,time,cell_id} = prpb_resp;
            end
        end
        cell_cnt{site,layer} = size(proj_meta(site).rd(layer,1).act,1);
    end
end

for site = 1:9
    for layer = 1:4
        if cell_cnt{site,layer} > 0
            for cell_id = 1:cell_cnt{site,layer}
                cur_fb_cell_act = vertcat(rps_resp_by_cell{site,layer,:,cell_id});
                if size(cur_fb_cell_act,1) == 1
                    cur_fb_cell_act = cur_fb_cell_act-mean(cur_fb_cell_act(1:7));
                    mean_rps_resp_by_cell{site,layer,cell_id} = cur_fb_cell_act;
                end
                if size(cur_fb_cell_act,1) > 1
                    mean_cur_fb_cell_act = mean(cur_fb_cell_act);
                    mean_cur_fb_cell_act = mean_cur_fb_cell_act-mean(mean_cur_fb_cell_act(1:7));
                    mean_rps_resp_by_cell{site,layer,cell_id} = mean_cur_fb_cell_act;
                end
            end
        end
    end
end

for site = 1:9
    for layer = 1:4
        if cell_cnt{site,layer} > 0
            for cell_id = 1:cell_cnt{site,layer}
                cur_pb_cell_act = vertcat(rpb_resp_by_cell{site,layer,:,cell_id});
                if size(cur_pb_cell_act,1) == 1
                    cur_pb_cell_act = cur_pb_cell_act-mean(cur_pb_cell_act(1:7));
                    mean_rpb_resp_by_cell{site,layer,cell_id} = cur_pb_cell_act;
                end
                if size(cur_pb_cell_act,1) > 1
                    mean_cur_pb_cell_act = mean(cur_pb_cell_act);
                    mean_cur_pb_cell_act = mean_cur_pb_cell_act-mean(mean_cur_pb_cell_act(1:7));
                    mean_rpb_resp_by_cell{site,layer,cell_id} = mean_cur_pb_cell_act;
                end
            end
        end
    end
end

resp_fb_act = vertcat(mean_rps_resp_by_cell{act_site,:,:});
resp_fb_pas = vertcat(mean_rps_resp_by_cell{pas_site,:,:});

mean_resp_fb_act = mean(resp_fb_act);
mean_resp_fb_pas = mean(resp_fb_pas);

sem_resp_fb_act = std(resp_fb_act)/sqrt(size(resp_fb_act,1));
sem_resp_fb_pas = std(resp_fb_pas)/sqrt(size(resp_fb_pas,1));

resp_pb_act = vertcat(mean_rpb_resp_by_cell{act_site,:,:});
resp_pb_pas = vertcat(mean_rpb_resp_by_cell{pas_site,:,:});

mean_resp_pb_act = mean(resp_pb_act);
mean_resp_pb_pas = mean(resp_pb_pas);

sem_resp_pb_act = std(resp_pb_act)/sqrt(size(resp_pb_act,1));
sem_resp_pb_pas = std(resp_pb_pas)/sqrt(size(resp_pb_pas,1));

figure;hold on
h1 = area([-7/15:1/15:75/15 75/15:-1/15:-7/15],[mean_resp_fb_act+sem_resp_fb_act fliplr(mean_resp_fb_act-sem_resp_fb_act)]*100);
h2 = area([-7/15:1/15:75/15 75/15:-1/15:-7/15],[mean_resp_fb_pas+sem_resp_fb_pas fliplr(mean_resp_fb_pas-sem_resp_fb_pas)]*100);
h3 = area([-7/15:1/15:75/15 75/15:-1/15:-7/15],[mean_resp_pb_act+sem_resp_pb_act fliplr(mean_resp_pb_act-sem_resp_pb_act)]*100);
h4 = area([-7/15:1/15:75/15 75/15:-1/15:-7/15],[mean_resp_pb_pas+sem_resp_pb_pas fliplr(mean_resp_pb_pas-sem_resp_pb_pas)]*100);

plot(-7/15:1/15:75/15,mean_resp_fb_act*100)
plot(-7/15:1/15:75/15,mean_resp_fb_pas*100,'r')

plot(-7/15:1/15:75/15,mean_resp_pb_act*100,'--')
plot(-7/15:1/15:75/15,mean_resp_pb_pas*100,'r--')

xlim([-0.5 5])
% ylim([-2 5])

set(h1,'FaceColor',[.5 .5 .5])
set(h2,'FaceColor',[.5 .5 .5])
set(h3,'FaceColor',[.5 .5 .5])
set(h4,'FaceColor',[.5 .5 .5])
set(h1,'LineStyle','none')
set(h2,'LineStyle','none')
set(h3,'LineStyle','none')
set(h4,'LineStyle','none')